Project Summary

This project aimed to look at the spatial variability of Symbiodinium clades C and D in the Kane’ohe Bay, O’ahu, Hawai’i population of Montipora capitata. We investigated the distribution of symbionts at scales ranging from location in the bay to location on an individual reef. We also looked at differences among reef types (fringing vs. patch), colony color morph (brown vs. orange) and depth. Heterogeneous mixtures of symbiont clades were considered in the analysis for spatial patterns. By investigating spatial variability of Symbiodinium, we furthered the understanding of stress-response potential in Kane’ohe Bay.

Library Packages

setwd("~/Symcap")
library(data.table)
library(lsmeans)
library(devtools)
library(plyr)
library(reshape2)
library(popbio)
library(RgoogleMaps)
library(plotrix)
library(zoo)
library(rgdal)
library(car)
library(scales)
library(png)
library(pixmap)
library(ecodist)
library(cluster)
library(fpc)
library(clustsig)
library(mapplots)
library(foreign)
library(nnet)
library(ggplot2)
library(mlogit)

Import/Merge Field and qPCR Data

Coral_Data <- read.csv("Coral_Collection.csv")
Coral_Data$Depth..m. <- as.numeric(as.character(Coral_Data$Depth..m.))
source_url("https://raw.githubusercontent.com/jrcunning/steponeR/master/steponeR.R")
Mcap.plates <- list.files(path = "qPCR_data", pattern = "txt$", full.names = T)
Mcap <- steponeR(files = Mcap.plates, delim="\t",
                 target.ratios=c("C.D"), 
                 fluor.norm = list(C=2.26827, D=0), 
                 copy.number=list(C=33, D=3))
Mcap <- Mcap$result
Mcap <- Mcap[grep("+", Mcap$Sample.Name, fixed=T, invert = T), ]
Mcap <- Mcap[grep("NTC", Mcap$Sample.Name, fixed = T, invert = T), ]
Mcap <- Mcap[grep("PCT", Mcap$Sample.Name, fixed = T, invert = T), ]
colnames(Mcap)[which(colnames(Mcap)=="Sample.Name")] <- "Colony"
Mcap$fail <- ifelse(Mcap$C.reps < 2 & Mcap$D.reps < 2, TRUE, FALSE)
fails <- Mcap[Mcap$fail==TRUE, ]
Mcap <- Mcap[which(Mcap$fail==FALSE),]
Mcap$C.D[which(Mcap$C.reps<2)] <- -Inf
Mcap$C.D[which(Mcap$D.reps<2)] <- Inf
Mcap <- Mcap[with(Mcap, order(Colony)), ]
Mcap$propC <- Mcap$C.D / (Mcap$C.D+1)
Mcap$propD <- 1-Mcap$propC
Mcap$propD[which(Mcap$C.D==-Inf)] <-1
Mcap$propC[which(Mcap$C.D==-Inf)] <-0
Mcap$propD[which(Mcap$C.D==Inf)] <-0
Mcap$propC[which(Mcap$C.D==Inf)] <-1
Mcap$Dom <- ifelse(Mcap$propC>Mcap$propD, "C", "D")
Symcap<-merge(Coral_Data, Mcap, by="Colony", all=T)
Symcap <- Symcap[with(Symcap, order(Colony)), ]
Symcap$Mix <- factor(ifelse(Symcap$propC>Symcap$propD, ifelse(Symcap$propD!=0, "CD", "C"), ifelse(Symcap$propD>Symcap$propC, ifelse(Symcap$propC!=0, "DC", "D"), NA)), levels = c("C", "CD", "DC", "D"))
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")

Change Bay Area for 18, 26 and F7-18

Symcap$Bay.Area[which(Symcap$Reef.ID=="26")] <- "Central"
Symcap$Bay.Area[which(Symcap$Reef.ID=="18")] <- "Southern"
Symcap$Bay.Area[which(Symcap$Reef.ID=="F7-18")] <- "Southern"

Adjust Depth by Mean Sea Level

To account for the influence of tides, depth was adjusted according to the difference in sea level from the mean sea level on each collection day at 6-minute intervals. Mean sea level was obtained from NOAA tide tables for Moku o Lo’e.

JuneTide=read.csv("Station_1612480_tide_ht_20160601-20160630.csv")
JulyTide=read.csv("Station_1612480_tide_ht_20160701-20160731.csv")
AugustTide=read.csv("Station_1612480_tide_ht_20160801-20160812.csv")
Tide<-rbind(JuneTide, JulyTide, AugustTide)
Tide$Time <- as.POSIXct(Tide$TimeUTC, format="%Y-%m-%d %H:%M:%S", tz="UTC")
attributes(Tide$Time)$tzone <- "Pacific/Honolulu"
Symcap$Time2 <- as.POSIXct(paste(as.character(Symcap$Date), as.character(Symcap$Time)),
                                format="%m/%d/%y %H:%M", tz="Pacific/Honolulu")
Symcap$Time=Symcap$Time2

# Add estimated times for missing values
Symcap$Time[170] <- as.POSIXct("2016-06-14 12:07:00")
Symcap$Time[177] <- as.POSIXct("2016-06-14 12:20:00")
Symcap$Time[178] <- as.POSIXct("2016-06-14 12:22:00")
Symcap$Time[180] <- as.POSIXct("2016-06-14 13:08:00")
Symcap$Time[187] <- as.POSIXct("2016-06-14 12:42:00")
Symcap$Time[188] <- as.POSIXct("2016-06-14 12:44:00")
Symcap$Time[206] <- as.POSIXct("2016-06-16 13:10:00")
Symcap$Time[208] <- as.POSIXct("2016-06-16 13:24:00")
Symcap$Time[211] <- as.POSIXct("2016-06-16 12:37:00")
Symcap$Time[218] <- as.POSIXct("2016-06-16 12:27:00")
Symcap$Time[448] <- as.POSIXct("2016-07-16 13:32:00")

Round6 <- function (times)  {
  x <- as.POSIXlt(times)
  mins <- x$min
  mult <- mins %/% 6
  remain <- mins %% 6
  if(remain > 3L) {
    mult <- mult + 1
  } else {
    x$min <- 6 * mult
  }
  x <- as.POSIXct(x)
  x
}

Symcap$Time.r <- Round6(Symcap$Time)
Tide$Time.r <- Tide$Time
merged<-merge(Symcap, Tide, by="Time.r", all.x=T)
merged$newDepth <- merged$Depth..m.- merged$TideHT

Coral Collection Map

#KB <- c(21.46087401, -157.809907) 
#KBMap <- GetMap(center = KB, zoom = 13, maptype = "satellite", SCALE = 2, GRAYSCALE = FALSE)
#save(KBMap, file = "KBMap.Rdata")
load("KBMap.Rdata")
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
XY <- subset(XY, Reef.ID!="37")
par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude, col=153, pch=21, bg="#FF7F50", lwd=2, cex = 1.2)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)
par(new=T, mar=c(14,21,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83") 
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83")) 
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()

Chi-Squared

Significant Differences

Dominant Symbiont by Reef Area

Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Dom, Symcap$Reef.Area)
chisq.test(results)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 136.26, df = 1, p-value < 2.2e-16
prop.table(results, margin = 2)
##    
##         Slope       Top
##   C 0.7767857 0.3294574
##   D 0.2232143 0.6705426
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Symbmiont Proportion")
legend("topright", legend=c("C", "D"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)

Dominant Symbiont by Color Morph

results=table(Symcap$Dom, Symcap$Color.Morph)
chisq.test(results)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 164.96, df = 1, p-value < 2.2e-16
prop.table(results, margin = 2)
##    
##         Brown    Orange
##   C 0.8896321 0.4103194
##   D 0.1103679 0.5896806
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Color Morph", ylab = "Symbiont Proportion")
legend("topright", legend=c("C", "D"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)

Symbiont Community Composition by Color Morph

results=table(Symcap$Mix, Symcap$Color.Morph)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 167.44, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##     
##            Brown      Orange
##   C  0.762541806 0.361179361
##   CD 0.127090301 0.049140049
##   DC 0.107023411 0.570024570
##   D  0.003344482 0.019656020
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Color Morph", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)

Symbiont Community Composition by Reef Area

Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Mix, Symcap$Reef.Area)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 138.97, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##     
##            Slope         Top
##   C  0.678571429 0.275193798
##   CD 0.098214286 0.054263566
##   DC 0.214285714 0.651162791
##   D  0.008928571 0.019379845
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Reef Area", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)

Symbiont Community Composition by Dominant Symbiont

results=table(Symcap$Mix, Symcap$Dom)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 707, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##     
##               C          D
##   C  0.86635945 0.00000000
##   CD 0.13364055 0.00000000
##   DC 0.00000000 0.96703297
##   D  0.00000000 0.03296703
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray 10", "gray 85", "gray 40", "gray100"), xlab = "Dominant Symbiont", ylab = "Symbiont Mixture Proportion")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray85", "gray40", "gray100"), inset = c(-.2, 0), xpd = NA)

When looking at D-only colonies, the Ct values are on the higher end of the spectrum (35 or greater) indicating poor amplification and the potential for C amplification being pushed back in the cycle. This is supported by the fact that 5 of the 9 D-only colonies had C present in 1 qPCR replicate.

Proportion of D When Present in Mixture

propD <- merged$propD[which(merged$propD > 0 & merged$propD < 1)]
hist(propD, xlab = "Proportion of Clade D", ylab = "Number of Samples", main = "", col = "gray75")

propDHist <- subset(merged, propD > 0 & propD < 1)
propDHist$propD <- cut(propDHist$propD, breaks = 5)
DCol=table(propDHist$Color.Morph, propDHist$propD)
z <- with(subset(merged, propD==0), table(Color.Morph))
o <- with(subset(merged, propD==1), table(Color.Morph))
DCol <- cbind(z, DCol, o)
par(mar=c(2, 4, 2, 0))
bars <- barplot(DCol, xlab = "Clade D Proportion", ylab = "Number of Samples", 
        main = "", col = c(alpha("sienna", 0.55), alpha("orange", 0.55)), 
        axisnames = FALSE, space = c(0,0.3,0,0,0,0,0.3))
axis(side = 1, at = c(1.3, 2.3, 3.3, 4.3, 5.3, 6.3), labels = c(">0.0", 0.2, 0.4, 0.6, 0.8, "<1.0"))
axis(side = 1, at = c(0.5, 7.2), labels = c("All C", "All D"), tick = FALSE)

propDHist <- subset(merged, propD > 0 & propD < 1)
propDHist$propD <- cut(propDHist$propD, breaks = 5)
DCol=table(propDHist$Color.Morph, propDHist$propD)
z <- with(subset(merged, propD==0), table(Color.Morph))
o <- with(subset(merged, propD==1), table(Color.Morph))
DCol <- cbind(z, DCol, o)
par(mar=c(2, 4, 2, 0))
bars <- barplot(DCol, xlab = "Clade D Proportion", ylab = "Number of Samples", 
        main = "", col = c(alpha("sienna", 0.75), alpha("orange", 0.75)), 
        axisnames = FALSE, space = c(0,0.3,0,0,0,0,0.3))
axis(side = 1, at = c(1.3, 2.3, 3.3, 4.3, 5.3, 6.3), labels = c(">0.0", 0.2, 0.4, 0.6, 0.8, "<1.0"))
axis(side = 1, at = c(0.5, 7.2), labels = c("All C", "All D"), tick = FALSE)

Color Morph by Reef Area

Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Color.Morph, Symcap$Reef.Area)
chisq.test(results)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 81.109, df = 1, p-value < 2.2e-16
prop.table(results, margin = 2)
##         
##              Slope       Top
##   Brown  0.5523385 0.2015504
##   Orange 0.4476615 0.7984496
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Color Morph Proportion")
legend("topright", legend=c("Brown", "Orange"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)

Color Morph by Bay Area

results=table(merged$Color.Morph, merged$Bay.Area)
prop.table(results, margin = 2)
##         
##            Central  Northern  Southern
##   Brown  0.4855769 0.3641026 0.4210526
##   Orange 0.5144231 0.6358974 0.5789474
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 6.1032, df = 2, p-value = 0.04728

Color Morph by Reef ID

results=table(Symcap$Color.Morph, Symcap$Reef.ID)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 50.691, df = 24, p-value = 0.001156
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propCol=table(Symcap$Color.Morph, Symcap$Reef.ID)
propCol=prop.table(propCol, margin = 2)
propCol <- as.data.frame.matrix(propCol)
props <- data.frame(t(propCol))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)

rownames(XY) <- XY$Reef.ID
XY <- XY[, -1]
XY <- na.omit(XY)
apply(XY, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["Orange"], reef["Brown"]), radius = 7, col = c("orange", "sienna"))
})

Mantel Test for Spatial Autocorrelation

Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propCol=table(Symcap$Color.Morph, Symcap$Reef.ID)
propCol=prop.table(propCol, margin = 2)
propCol <- as.data.frame.matrix(propCol)
props <- data.frame(t(propCol))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
reef.dists <- dist(cbind(XY$Longitude, XY$Latitude))
col.dists <- dist(XY$Brown)
set.seed(12456)
mantel(col.dists~reef.dists)
##     mantelr       pval1       pval2       pval3   llim.2.5%  ulim.97.5% 
## -0.09132955  0.87800000  0.12300000  0.24500000 -0.14828356 -0.02809549

Color Morph by Bay Area Adjusted for Depth

merged$Color <- ifelse(merged$Color.Morph=="Orange", 1, 0)
dat <- aggregate(data.frame(prop=merged$Color), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) 
mod <- glm(Color ~ newDepth + Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "Raw", "O.Adj")

Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propCol=table(Symcap$Color.Morph, Symcap$Reef.ID)
propCol=prop.table(propCol, margin = 2)
propCol <- as.data.frame.matrix(propCol)
props <- data.frame(t(propCol))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$B.Adj <- 1-XY2$O.Adj

PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)

rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["O.Adj"], reef["B.Adj"]), radius = 7, 
               col = c("orange", "sienna"))
})

Mantel Test for Spatial Autocorrelation

reef.dists <- dist(cbind(XY$Longitude, XY$Latitude))
col.dists <- dist(XY2$O.Adj)
set.seed(12456)
mantel(col.dists~reef.dists)
##     mantelr       pval1       pval2       pval3   llim.2.5%  ulim.97.5% 
## -0.07857358  0.84000000  0.16100000  0.31400000 -0.12666742 -0.01919027

Symbiont Community Composition by Bay Area

results=table(Symcap$Mix, Symcap$Bay.Area)
prop.table(results, margin = 2)
##     
##         Central   Northern   Southern
##   C  0.50961538 0.61538462 0.49174917
##   CD 0.08173077 0.05641026 0.09900990
##   DC 0.39423077 0.30256410 0.40594059
##   D  0.01442308 0.02564103 0.00330033
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 14.719, df = 6, p-value = 0.02256

Non-Significant Differences

Dominant Symbiont by Reef Type

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 0.77593, df = 1, p-value = 0.3784

Dominant Symbiont by Bay Area

results=table(Symcap$Dom, Symcap$Bay.Area)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 3.8852, df = 2, p-value = 0.1433

Dominant Symbiont by Reef ID

results=table(Symcap$Dom, Symcap$Reef.ID)
chisq.test(results)
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propDom=table(Symcap$Dom, Symcap$Reef.ID)
propDom=prop.table(propDom, margin = 2)
propDom <- as.data.frame.matrix(propDom)
props <- data.frame(t(propDom))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)

rownames(XY) <- XY$Reef.ID
XY <- XY[, -1]
XY <- na.omit(XY)
apply(XY, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["C"], reef["D"]), radius = 7, col = c("blue", "red"))
})

Mantel Test for Spatial Autocorrelation

Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propDom=table(Symcap$Dom, Symcap$Reef.ID)
propDom=prop.table(propDom, margin = 2)
propDom <- as.data.frame.matrix(propDom)
props <- data.frame(t(propDom))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
reef.dists <- dist(cbind(XY$Longitude, XY$Latitude))
dom.dists <- dist(XY$C)
set.seed(12456)
mantel(dom.dists~reef.dists)
##    mantelr      pval1      pval2      pval3  llim.2.5% ulim.97.5% 
## 0.15974663 0.04000000 0.96100000 0.04500000 0.08875833 0.24328773

pval1 indicates that there is a positive correlation between distance and difference among reefs in terms of dominant symbiont clade. Reefs that are closer in proximity to each other are more similar in terms of symbiont dominance. A positive mantelr value indicates a positive correlation.

Clustering of Similar Reefs

doms <- hclust(dom.dists)
plot(doms, labels = XY$Reef.ID)

a <- hclust(reef.dists*dom.dists)
plot(a, labels = XY$Reef.ID)

domsk <- pamk(data = dom.dists, critout = TRUE)
domsk

km <- kmeans(dom.dists, centers = 5)
km
df <- cbind(XY, km$cluster)
km$cluster

PlotOnStaticMap(KBMap, df$Latitude, df$Longitude, col="red", pch=as.character(km$cluster), lwd = 2)

Dominant Symbiont Plots

Raw Values
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
XY <- subset(XY, Reef.ID!="37")

merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")

XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c

PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)

rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["raw.c"], reef["raw.d"]), radius = 7, 
               col = c("red", "blue"))
})

Depth Adjustment
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")

XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c

PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)

rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["c.adj"], reef["d.adj"]), radius = 7, 
               col = c("red", "blue"))
})

Depth Interaction Adjustment
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")

XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c

PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)

rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["c.adj.int"], reef["d.adj.int"]), radius = 7, 
               col = c("red", "blue"))
})

Depth Adjustment for Color Morph and Dominant Symbiont

A multinomial logistic regression was performed on the interaction of color morph and dominant symbiont to discount the influence of depth on spatial distribution.

merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), 
                 FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")

XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c

manteltable = table(merged$Dom, merged$Color.Morph, merged$Reef.ID)
nc <- aggregate(interaction(merged$Color.Morph, merged$Dom), 
                by=list(merged$Reef.ID), FUN=table)
nc <- data.frame(Reef.ID=as.character(nc[,1]), prop.table(nc[,2], margin=1))
XY3 <- merge(XY2, nc, by="Reef.ID", all=T)

merged$ColDom <- interaction(merged$Color.Morph, merged$Dom)
model <- multinom(ColDom~newDepth+Reef.ID, merged)
means <- lsmeans(model, specs = "Reef.ID")
means
pp <- fitted(model)

newdat <- data.frame(Reef.ID = levels(merged$Reef.ID),
                   newDepth = mean(merged$newDepth, na.rm=T))
pred <- predict(model, newdata = newdat, "probs")
new <- data.frame(Reef.ID=as.character(newdat[,1]), pred)
XY4 <- merge(XY3, new, by="Reef.ID", all=T)

PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)

XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["Brown.C.x"], reef["Orange.C.x"], reef["Brown.D.x"], 
                   reef["Orange.D.x"]), radius = 7, 
   col = c("turquoise4", "steelblue1", "tomato", "coral"))
})

legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("turquoise4", "steelblue1", "tomato", "coral"))

par(new=T, mar=c(15,22,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83") 
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83")) 
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()

Raw proportion values for Color Morph and Dominant Symbiont at each reef not discounting the effect of depth.

Mantel Test for Spatial Autocorrelation
reef.dists <- dist(cbind(XY4$Longitude, XY4$Latitude))
dom.dists <- bcdist(cbind(XY4$Brown.C.x, XY4$Orange.C.x, XY4$Brown.D.x, XY4$Orange.D.x))
set.seed(12456)
mantel(dom.dists~reef.dists)
##    mantelr      pval1      pval2      pval3  llim.2.5% ulim.97.5% 
## 0.09077084 0.12000000 0.88100000 0.24000000 0.01812853 0.14910712
PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)

XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["Brown.C.y"], reef["Orange.C.y"], reef["Brown.D.y"], 
                   reef["Orange.D.y"]), radius = 7, 
               col = c("turquoise4", "steelblue1", "tomato", "coral"))
})

legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("turquoise4", "steelblue1", "tomato", "coral"))

par(new=T, mar=c(15,22,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83") 
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83")) 
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()

Adjusted values of Dominant Symbiont and Color Morph adjusted for the influence of depth.

Mantel Test for Spatial Autocorrelation
reef.dists <- dist(cbind(XY4$Longitude, XY4$Latitude))
dom.dists <- bcdist(cbind(XY4$Brown.C.y, XY4$Orange.C.y, XY4$Brown.D.y, XY4$Orange.D.y))
set.seed(12456)
mantel(dom.dists~reef.dists)
##      mantelr        pval1        pval2        pval3    llim.2.5% 
##  0.048985367  0.229000000  0.772000000  0.493000000 -0.008775995 
##   ulim.97.5% 
##  0.110127704
Color and Dominant Symbiont by Bay Area

When considering the effect of bay area on the interaction of color morph and dominant symbiont, the proportion of Brown colonies dominated by D increases as you move from the north to the south of the bay. The proportion increases almost 6x, yet the small number of colonies makes this non-compelling.

Type=table(merged$ColDom, merged$Bay.Area)
prop.table(Type, margin = 2)
##           
##               Central   Northern   Southern
##   Brown.C  0.42307692 0.35384615 0.35973597
##   Orange.C 0.16826923 0.31794872 0.23102310
##   Brown.D  0.06250000 0.01025641 0.05940594
##   Orange.D 0.34615385 0.31794872 0.34983498
chisq.test(Type)
## 
##  Pearson's Chi-squared test
## 
## data:  Type
## X-squared = 19.376, df = 6, p-value = 0.003573

Color Morph by Reef Type

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 1.1302, df = 1, p-value = 0.2877

Symbiont Community Composition by Reef Type

## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 2.8371, df = 3, p-value = 0.4174

Symbiont Community Composition by Reef ID

## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 87.127, df = 72, p-value = 0.1081

Logistic Regression

Significant Effect

Dominant Symbiont by Depth

merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
Dom1 <- subset(merged, !is.na(newDepth) & !is.na(Dominant))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Dominant
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                       705     942.15              
## newDepth  1   129.01       704     813.13 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
logi.hist.plot(Dom1$newDepth, Dom1$Dominant, boxp = FALSE, type = "hist", col="gray", xlabel = "Depth (m)", ylabel = "", ylabel2 = "")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "C                             D", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of clade C Symbiont", line = 3, cex = 1)

merged$DepthInt <- cut(merged$newDepth, breaks = 0:13)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results
##    
##     (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
##   0   149    33    27    23    14     4     3     2     2      1       0
##   1    73    52    83    92    46    27    15    12    16      8       5
##    
##     (11,12] (12,13]
##   0       1       0
##   1       0       1
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.2, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Probability of D-Dominance")

Bars indicate relative frequency of clade C vs. D dominance at 1m depth intervals when pooling all samples collected.

Mixture by Depth

merged$Mixture <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 1, 0)
merged$Mixture2 <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 0, 1)
results=glm(Mixture~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Mixture
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                       705     973.27              
## newDepth  1   93.723       704     879.55 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
results=table(merged$Mixture2, merged$DepthInt)
results
##    
##     (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
##   0   154    44    39    28    18     8     4     2     5      2       1
##   1    68    41    71    87    42    23    14    12    13      7       4
##    
##     (11,12] (12,13]
##   0       1       0
##   1       0       1
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Mix", "Non-Mix"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.23, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Mixture~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Mixture Proportion")

Bars indicate relative frequency of Mixture vs. Non-Mixtures at 1m depth intervals when pooling all samples collected.

Color Morph by Depth

merged$Color <- ifelse(merged$Color.Morph=="Orange", 1, 0)
results=glm(Color~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Color
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                       706     963.85              
## newDepth  1   78.431       705     885.42 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Color <- subset(merged, !is.na(newDepth) & !is.na(Color))
logi.hist.plot(independ = Color$newDepth, depend = Color$Color, type = "hist", boxp = FALSE, ylabel = "", col="gray", ylabel2 = "", xlabel = "Depth (m)")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "Brown                       Orange", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)

merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results
##    
##     (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
##   0   171    42    57    63    30    11     9     2     3      1       0
##   1    51    43    53    52    30    20     9    12    16      8       5
##    
##     (11,12] (12,13]
##   0       1       0
##   1       0       1
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(-.22, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Probability of Orange-Dominance")

Bars indicate relative frequency of Brown vs. Orange color morph dominance at 1m depth intervals when pooling all samples collected.

merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results
##    
##     (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
##   0   171    42    57    63    30    11     9     2     3      1       0
##   1    51    43    53    52    30    20     9    12    16      8       5
##    
##     (11,12] (12,13]
##   0       1       0
##   1       0       1
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.45), alpha("sienna", 0.45)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.45), alpha("orange", 0.45)), inset = c(-.22, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Probability of Orange-Dominance")

Two-Way ANOVA

Significant Interactive Effects

Dominant Symbiont by Color Morph and Reef Area

merged$Reef.Area <- ifelse(merged$Reef.Area!="Top", yes = "Slope", no = "Top")
table(merged$Dom, merged$Color.Morph, merged$Reef.Area)
## , ,  = Slope
## 
##    
##     Brown Orange
##   C   226    122
##   D    21     79
## 
## , ,  = Top
## 
##    
##     Brown Orange
##   C    40     45
##   D    12    161
model=aov(Dominant~Color.Morph*Reef.Area, data = merged)
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: Dominant
##                        Sum Sq  Df F value    Pr(>F)    
## Color.Morph            21.329   1 134.208 < 2.2e-16 ***
## Reef.Area              14.489   1  91.168 < 2.2e-16 ***
## Color.Morph:Reef.Area   1.780   1  11.201 0.0008612 ***
## Residuals             111.565 702                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model=aov(Dominant~newDepth*Bay.Area, data = merged)
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: Dominant
##                    Sum Sq  Df  F value  Pr(>F)    
## newDepth           25.026   1 125.1913 < 2e-16 ***
## Bay.Area            0.954   2   2.3870 0.09265 .  
## newDepth:Bay.Area   1.555   2   3.8886 0.02092 *  
## Residuals         139.933 700                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dominant Symbiont by Color Morph and Depth

Because an interactive effect of reef area and color morph was observed and slope is indicative of a depth gradient, the interaction between depth and color morph was tested here.

model1=lm(Dominant~Color.Morph*newDepth, data = merged)
Anova(model1, type = 2)
## Anova Table (Type II tests)
## 
## Response: Dominant
##                       Sum Sq  Df F value    Pr(>F)    
## Color.Morph           24.385   1 152.958 < 2.2e-16 ***
## newDepth               9.778   1  61.334 1.782e-14 ***
## Color.Morph:newDepth   6.140   1  38.514 9.290e-10 ***
## Residuals            111.916 702                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Depth at which Orange switches from D to C Dominance

While brown was always C-dominated, orange was observed to switch from D to C dominance. The depth threshold at which this switch occurs was analyzed here.

threshdepth <- function(color) {
  df <- subset(merged, Color.Morph==color)
  plot(df$Dominant2~df$newDepth, xlab="Depth (m)", ylab = "Probability of Clade C Symbiont",
       main=color)
  abline(h = 0.5, lty=2)
  results=glm(Dominant2~newDepth, family = "binomial", data = df)
  pval <- data.frame(coef(summary(results)))$`Pr...z..`[2]
  mtext(side=3, text=pval)
  newdata <- list(newDepth=seq(0,11,0.01))
  fitted <- predict(results, newdata = newdata, type = "response")
  lines(fitted ~ seq(0,11,0.01))
  thresh <- ifelse(pval < 0.05,
                   newdata$newDepth[which(diff(sign(fitted - 0.5))!=0)], NA)
  return(thresh)
}
sapply(levels(merged$Color.Morph), FUN=threshdepth)

##  Brown Orange 
##     NA   2.75
df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,11,0.01), col="sienna", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Clade D Symbiont", line = 3, cex = 1)
legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(-.22, 0), xpd = NA)

Dominant Symbiont by Depth and Reef Type

model3=aov(Dominant2~newDepth*Reef.Type, data = merged)
Anova(model3, type = 2)
## Anova Table (Type II tests)
## 
## Response: Dominant2
##                     Sum Sq  Df  F value    Pr(>F)    
## newDepth            24.785   1 123.5859 < 2.2e-16 ***
## Reef.Type            0.014   1   0.0674  0.795285    
## newDepth:Reef.Type   1.643   1   8.1944  0.004327 ** 
## Residuals          140.785 702                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df <- subset(merged, Reef.Type=="Patch")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,11,0.01), col="brown1", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Clade D Symbiont", line = 3, cex = 1)
legend("topright", legend=c("Patch", "Fringe"), fill=c("dodgerblue3", "brown1"), inset = c(-.22, 0), xpd = NA)

Color Morph by Depth and Reef Type

model2=aov(Color~newDepth*Reef.Type, data = merged)
Anova(model2, type = 2)
## Anova Table (Type II tests)
## 
## Response: Color
##                     Sum Sq  Df F value  Pr(>F)    
## newDepth            18.186   1 83.6137 < 2e-16 ***
## Reef.Type            0.065   1  0.2993 0.58448    
## newDepth:Reef.Type   1.294   1  5.9504 0.01496 *  
## Residuals          152.900 703                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df <- subset(merged, Reef.Type=="Patch")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,11,0.01), col="brown1", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)
legend("topright", legend=c("Patch", "Fringe"), fill=c("dodgerblue3", "brown1"), inset = c(-.22, 0), xpd = NA)

Dominant Symbiont per Color Morph by Depth and Reef Type

Because depth and reef type have an interactive effect on both color morph and dominant symbiont clade, the dominant symbiont per color morph at each reef type was visualized.

df <- subset(merged, Reef.Type=="Patch")
dfo <- subset(df, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = dfo)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)

df <- subset(merged, Reef.Type=="Patch")
dfb <- subset(df, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = dfb)
newdata <- list(newDepth=seq(0,11,0.01))
par(new = T, mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="black", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)

df <- subset(merged, Reef.Type=="Fringe")
dfo <- subset(df, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = dfo)
newdata <- list(newDepth=seq(0,11,0.01))
par(new = T, mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="red", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)

df <- subset(merged, Reef.Type=="Fringe")
dfb <- subset(df, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = dfb)
newdata <- list(newDepth=seq(0,11,0.01))
par(new = T, mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
legend("topright", legend=c("PO", "PB", "FO", "FB"), fill=c("dodgerblue3", "black", "red", "orange"), inset = c(-.2, 0), xpd = NA)
mtext(side = 2, text = "Probability of Clade D Symbiont", line = 3, cex = 1)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)

Dominant Symbiont by Depth and Reef ID

model2=aov(Dominant~newDepth*Reef.ID, data = merged)
Anova(model2, type = 2)
## Anova Table (Type II tests)
## 
## Response: Dominant
##                   Sum Sq  Df  F value    Pr(>F)    
## newDepth          23.693   1 122.6080 < 2.2e-16 ***
## Reef.ID            6.905  24   1.4888  0.063088 .  
## newDepth:Reef.ID   8.772  24   1.8915  0.006452 ** 
## Residuals        126.765 656                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
merged$Dominant2 <- ifelse(merged$Dom=="C", 0, 1)
par(mfrow=c(5,5))

domreef <- function(id) {
  df <- subset(merged, Reef.ID==id)
  results=glm(Dominant2~newDepth, family = "binomial", data = df)
  newdata <- list(newDepth=seq(0,11,0.01))
  par(mar=c(1, 1, 3, 1))
  fitted <- predict(results, newdata = newdata, type = "response")
  plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", 
       col="dodgerblue3", lwd=3, xlab="", ylab="")
  mtext(side = 3, text = id)
  abline(h=0.5, lty=2)
}

sapply(levels(merged$Reef.ID), FUN=domreef)

On the y-axis, a value of 1 indicates D-dominance and a value of 0 indicates C-dominance.

Mixture by Depth and Reef ID

merged$Mixture <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 1, 0)
model4=aov(Mixture~newDepth*Reef.ID, data = merged)
Anova(model4, type = 2)
## Anova Table (Type II tests)
## 
## Response: Mixture
##                   Sum Sq  Df F value    Pr(>F)    
## newDepth          18.913   1 90.9371 < 2.2e-16 ***
## Reef.ID            9.987  24  2.0008  0.003245 ** 
## newDepth:Reef.ID   8.020  24  1.6068  0.033901 *  
## Residuals        136.434 656                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow=c(5,5))

mixreef <- function(id) {
  df <- subset(merged, Reef.ID==id)
  results=glm(Mixture~newDepth, family = "binomial", data = df)
  newdata <- list(newDepth=seq(0,11,0.01))
  par(mar=c(1, 1, 3, 1))
  fitted <- predict(results, newdata = newdata, type = "response")
  plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", 
       col="dodgerblue3", lwd=3, xlab="", ylab="")
  mtext(side = 3, text = id)
  abline(h=0.5, lty=2)
}

sapply(levels(merged$Reef.ID), FUN=mixreef)

On the y-axis, a value of 1 indicates a mixture and a value of 0 indicates single symbiont clade.

Figures

Orange and Brown Color Morph

Symbiont Community Composition by Dominant Symbiont

results=table(Symcap$Mix, Symcap$Dom)
chisq.test(results)
prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c(alpha("blue", 0.75), alpha("blue", 0.25), alpha("red", 0.25), alpha("red", 0.75)), xlab = "Dominant Symbiont", ylab = "Symbiont Mixture Proportion")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c(alpha("blue", 0.75), alpha("blue", 0.25), alpha("red", 0.25), alpha("red", 0.75)), inset = c(-.2, 0), xpd = NA)

When looking at D-only colonies, the Ct values are on the higher end of the spectrum (35 or greater) indicating poor amplification and the potential for C amplification being pushed back in the cycle. This is supported by the fact that 5 of the 9 D-only colonies had C present in 1 qPCR replicate.

Proportion of Clade D when Present

propD <- merged$propD[which(merged$propD > 0 & merged$propD < 1)]

propDHist <- subset(merged, propD > 0 & propD < 1)
propDHist$propD <- cut(propDHist$propD, breaks = 5)
DCol=table(propDHist$Color.Morph, propDHist$propD)
z <- with(subset(merged, propD==0), table(Color.Morph))
o <- with(subset(merged, propD==1), table(Color.Morph))
DCol <- cbind(z, DCol, o)
par(mar=c(4, 4, 2, 0))
bars <- barplot(DCol, xlab = "Clade D Proportion", ylab = "Number of Samples", main = "", col = c(alpha("sienna", 0.55), alpha("orange", 0.55)), axisnames = FALSE, space = c(0,0.3,0,0,0,0,0.3))
axis(side = 1, at = c(1.3, 2.3, 3.3, 4.3, 5.3, 6.3), labels = c(">0%", "20%", "40%", "60%", "80%", "<100%"))
axis(side = 1, at = c(0.5, 7.2), labels = c("All C", "All D"), tick = FALSE)

Depth Influence on Dominant Symbiont and Color Morph

par(mfrow=c(3,1))

merged$DepthInt <- cut(merged$newDepth, breaks = 0:13)
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(2, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(2.1, 4, 2, 6))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="Probability of D-Dominance")

merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(3, 4, 1, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)), 
        xlab = "", ylab = "Probability of Orange-Dominance",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(3.1, 4, 1, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="")

df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 0, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="Depth (m)", ylab="Probabilty of D-Dominance", axisnames=FALSE, xaxs = "i", yaxs = "i")
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0, 11, 0.01), col="sienna", lwd=3)
legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(0, 0), xpd = NA)

#results=table(merged$Dominant2, merged$Color.Morph)
#chisq.test(results)
#prop.table(results, margin = 2)
#par(new=T, mar=c(10, 10, .5, 6.3))
#barplot(prop.table(results, margin = 2), col = c(alpha("red", 0.35), alpha("blue", 0.35)), xlab = "", ylab = "", yaxt = 'n')
#legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.35), alpha("red", 0.35)), inset = c(0, 0), xpd = NA)

Color Morph and Dominant Symbiont at each Reef

par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)

XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["Brown.C.y"], reef["Orange.C.y"], reef["Brown.D.y"], 
                   reef["Orange.D.y"]), radius = 7, 
               col = c("royalblue3", "cornflowerblue", "red", "rosybrown2"))
})

legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("royalblue3", "cornflowerblue", "red", "rosybrown2"))

par(new=T, mar=c(14,21,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83") 
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83")) 
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()

Adjusted values of Dominant Symbiont and Color Morph adjusted for the influence of depth.

Color Morph and Dominant Symbiont per Bay Area

Latitude=aggregate(Latitude~Bay.Area, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Bay.Area, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Bay.Area", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY

merged$ColDom <- interaction(merged$Color.Morph, merged$Dom)
model <- multinom(ColDom~newDepth+Bay.Area, merged)
means <- lsmeans(model, specs = "Bay.Area")
means
newdat <- data.frame(Bay.Area = levels(merged$Bay.Area),
                   newDepth = mean(merged$newDepth, na.rm=T))
pred <- predict(model, newdata = newdat, "probs")
new <- data.frame(Bay.Area=as.character(newdat[,1]), pred)

XY4 <- merge(XY, new, by="Bay.Area", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Bay.Area

par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)

XY4$Bay.Area <- as.numeric(XY4$Bay.Area)
apply(XY4, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["Brown.C"], reef["Orange.C"], reef["Brown.D"], 
                   reef["Orange.D"]), radius = 20, 
               col = c("blue", "cornflowerblue", "red", "lightpink"))
})
legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("blue", "cornflowerblue", "red", "lightpink"))

par(new=T, mar=c(14,21,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83") 
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83")) 
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()

par(mar=c(3,3,0,0))
plot(HI, xlim=c(-157.82, -157.82), ylim=c(21.42, 21.51), 
     lwd=2, col="gray") 
legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), 
       fill = c("turquoise4", "steelblue1", "tomato", "coral"), bg = "white")
add.pie(z = c(42, 23, 5, 30), x = -157.8109, y = 21.46222, labels = NA, radius = 0.005, 
        col = c("turquoise4", "steelblue1", "tomato", "coral"))
add.pie(z = c(36, 33, 2, 29), x = -157.8286, y = 21.47748, labels = NA, radius = 0.005,
        col= c("turquoise4", "steelblue1", "tomato", "coral"))
add.pie(z = c(40, 24, 8, 28), x = -157.7965, y = 21.44077, labels = NA, radius = 0.005,
        col= c("turquoise4", "steelblue1", "tomato", "coral"))
box()
par(new=T, mar=c(17,24,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83") 
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83")) 
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()